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Twin peaks kHz QPOs: mathematics of the 3 : 2 orbital resonance 
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Abstract 

Using the method of multiple scales, one can derive an analytic solution that describes the behaviour of 
weakly coupled, non-linear oscillations in nearly Keplerian discs around neutron stars or black holes close 
to the 3:2 orbital epicyclic resonance. The solution obtained agrees with the previous numerical simulation. 
Such result may be relevant to the kilohertz quasi-periodic variability in X-ray fluxes observed in many 
Galactic black hole and neutron star sources. With a particular choice of tunable parameters, the solution 
fits accurately the observational data for Sco X-1. 

Key words: General relativity - Accretion - X-rays: individual (Sco X-1) - QPOs 



> ■ 

^ ■ 

m ■ 
m ■ 

o : 
^ . 
o ■ 

Oh- 
I 

o ■ 



1. Introduction 

Many Galactic black hole and neutron star sources in 
low mass X-ray binaries show both chaotic and quasi- 
periodic variability in their observed X-ray fluxes. Some 
of the quasi-periodic oscillations (QPOs) are in the kHz 
range and often come in pairs (j^upp,^'down) of twin peaks in 
the Fourier power spectra (e.g., van der KHs, 2000). In all 
four black hole sources with twin peak kHz QPOs pairs. 



'upp 



-- 3/2 (McClintock & RemiUard, 2003). For 
neutron stars sources the ratio of twin-peak frequencies 
is mostly close to 3/2 too (Abramowicz, Bulik, Bursa, 
Kluzniak, 2003). Based on these and other properties 
of QPOs Kluzniak & Abramowicz (2000) concluded that 
twin peak kHz QPOs are due to a resonance in the ac- 
cretion disc oscillation modes and they noticed that the 
specific 3 : 2 ratio could be a consequence of strong gravity. 

According to the standard Shakura-Sunyaev accretion 
disc model, matter spirals down into the central black 
hole along stream lines that are located almost on the 
equatorial plane 9 = 9o = tt/2, and that locally differ only 
slightly from a family of concentric circles r = rp = const. 
The small deviations, Sr = r — rQ, Sd ~9 — 9o are governed, 
with accuracy to linear terms, by 



Sr + cij^ 5r ~ Sttr 



LuU9 = < 



(1) 



Here the dot denotes the time derivative. For purely 
Keplerian (free) motion 6ar = 0, Sag = and the above 
equations describe two uncoupled harmonic oscillators 
with the eigenfrequcncies cug, LUr- 

We shall start with an argument which appeals to phys- 
ical intuition and which shows that the resonance now 
discussed is a very natural, indeed necessary, consequence 
of strong gravity. Assume that in thin discs, random fluc- 
tuations have Sr 3> S9.^ Thus, Sr59 is a first order term 
in S9 and should be included in the first order equation 
for vertical oscillations (1). The equation now takes the 

^ In actual calculations this additional assumption is not made. 



form, 

S9 + Lj^[l + hSr]S9 = Sae, (2) 

where ft, is a known constant. The first order equation for 
Sr has the solution Sr = Aocos{ojrt). Inserting this in (2) 
together with Sag = 0, one arrives to the Mathieu equation 
{Aq is absorbed in h), 

S9 + uj^[l + hcos{tOrt)]5e = 0, (3) 

that describes the parametric resonance. From the theory 
of the Mathieu equation one knows that when 

uj,,. V,. 2 

UJ0 ve 

the parametric resonance is excited (Landau & Lifshitz, 
1976). The resonance is strongests for the smallest pos- 
sible value of n. Because near black holes and neutron 
stars Vr < vq (see Figure 1), the smallest possible value 
for resonance is n = 3, which means that 2vg = 3i'r- This 
is an example of a situation in which parametric resonance 
works in a thin accretion disc and it intuitively explains 
the observed 3:2 ratio (Kluzniak and Abramowicz, 2002), 
because, obviously, 

(5) 



— = — = ?i = l,2,3. 



(4) 



^upp — ^9 1 ^down 

Of course, in real discs neither Sr = AQCOsitUrt), nor 
Sag = exactly; moreover in realistic situations, for purely 
geodesic motion (Sag = Sar = 0) the system does not show 
increasing amplitudes (the higher terms prevent this from 
occurring) . 

However, one may expect that because these equations 
are approximately obeyed in thin discs, the parametric 
resonance will indeed be excited. Such a resonance was 
found in numerical simulations of oscillations in a nearly 
Keplerian accretion disc by Abramowicz et al. (2003). 

2. Nearly-geodesic motion 

The equations roughly outlined in the previous Section 
will now be derived. 
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Fig. 1. Epicyclic frequencies in Schwarzschild's metric: 
meridional (solid line) and radial one (dashed line). The ra- 
dius is in units of tq. The vertical lines indicates the radii 
{^3:2,r$:3,r4:2) where the ratio between the two frequencies is 
3:2,5:3 and 4 : 2 (from right to left). 

The equations of the geodesic motion of a unit mass 
test-particle in spherical coordinates can be written in the 
form {t = T /ro): 



(r) + 



dU,ffir,0) 



de 



+ 2rV(r)^(T) = 



(6) 



dr 



(r) +r;,f(T)^ = o 



0(r)=O , 

where Ueff = g** + Ig*"^ + I'^g'^'^ is the effective potential 
and E = —ut is the energy. 

Consider a geodesic on the equatorial plane (j'Oi^o = 
7r/2,0 = f2f) and perturb it slightly, keeping constant the 
angular momentum. 



'^{t') =ro + p{f) p«ro 
9{f) = ^ + z(f} z<<7r/2 



(7) 



Here p{f) and z(f ) denote small devations from the circu- 
lar orbit, and to first order they describe two uncoupled 
harmonic oscillators with epicyclic eigenfrequencies: 



1 d^a 



1 d^a 



eff 



The ratio of the two frequencies is 



> 1 



9rr dr /f,ro,7r/2 



(10) 



This does not happen in the Newtonian case where the 
radial and vertical frequencies are always equal to the or- 
bital frequency (Jl = GM/rp). When continuing to higher 
orders in the Taylor expansion, the two motions become 
coupled, and this is why the so-called parametric reso- 
nance can arise. The perturbed equations are nonlinear 
differential equations with constant coefficients which are 



the values of the different derivatives of the effective po- 
tential Uef / and the metric at equilibrium (see Appendix 
2): 

Sag = z{f)+ujgz{f) - f{p{f),z{f),ro,9o) (11) 
Sar = p{f) + uj'^^p{f) - g{p{f),z{f),ro,0Q) . 

In the case of Schwarzschild's metric, 
f{p,z,ro,9o) = ciizp + Cbzp+ (12) 

+ C2lP^Z + Clbpzp + C03Z^ 

g{p,z,ro,9o) = 6022^ + 620/9^ + 6222^ + 

+e30P^ + eize2pz^ + ei2pz^ + 

6r2P^ + eire2P^P ■ 

The third order terms in the Taylor expansion provide the 
damping which saturates the increasing amplitudes: this 
motivates stopping the expansion at this order of the per- 
turbation. 

In first approximation wc used Paczyhski & Wiita (1980) 
pseudo-Newtonian model. The perturbed equations in the 
model and in general relativity look the same, apart from 
those terms in g{p,z,rQ,9Q) which involve p^ that here are 
absent. 

On the geodesies 5ag — Sa^ — 0; however, in order to ex- 
plore the features of the resonance, a small additional 
isotropic force, parametrized by \a\ G [0,1], can be added 
(Abramowicz et al. 2003); this force at the moment has 
no physical meaning (it may be a coupling due for exam- 
ple to pressure or viscosity), it is just a mathematical tool. 
Then 



(13) 



Sag = -a(c2ip^ + 6032^)2 
5ar^ -a{e^2 + eize2p)z^ ■ 

When |a| = the previous equations are just the ap- 
proximation of the geodesies, while when |q;| increases, 
they describe slight deviations from the free motion. 



(8) 3. The method of multiple scales 



This method is a variation of the straightforward expan- 
sion which leads to a unformly valid approximate solution 
of systems of weakly nonlinear differential equations; the 
^j^dcrlying idea is to consider the expansion which repre- 
sents the approximate solution to be a function of mul- 
tiple independent variables, or scales, instead of a single 
variable. The new "time-like" independent variables are 
defined: Tk = e^t for k = 0,1,2... Expressing the solution 
as a function of more variables, treated as independent, is 
an artifice to remove the secular terms, which would make 
the solution to explode unphysically. 
By writing Tfc one makes a formal assumption of physi- 
cal slow variation explicit. Indeed define progressively 
longer time scales, which are not negligible when t is of 
the order of 1/6*^ or longer. The characteristic time scale 
of the orbital motion near a neutron star or a stellar mass 
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Fig. 2. The dashed line is the least-squares best-fit to the 
data points (the observed kHz QPOs frequencies in Sco 
X — 1); the thin solid line corresponds to a slope of 3 : 2 
(for reference) . The thick segment is the analytic approx- 
imation, in which the frequencies are scaled for comparison 
with observations. 

black hole is short (Tq ~ Ims): hence for e of the order of 
10^^ the dynamical effects on the scale T3 become impor- 
tant in the interval of Is. 

The conditions to eliminate the secular terms give a sim- 
pler system of differential equations: for a set of ini- 
tial conditions the solution to such a system is a sum 
of trigonometric functions. The same solution can be 
found by assuming a priori that the asymptotic solution is 
the sum of trigonometric functions in which the dominant 
terms have frequencies not far from the original eigcnfrc- 
qucncics of the uncoupled harmonic oscillators (method 
of Lindstcdt-Poincarc). The corrections to the cigenfrc- 
quencies depend on the amplitudes and on the constant 
coefficients Hence in a Fourier transform of the geodesies 
perturbed in the vicinity of r^-^ =27/5, two peaks centered 
near 3 : 2 would be found, the position of which and the 
ratio depend on the perturbation (a) and on the initial 
amplitudes. 

This method was used to find the approximate solution 
to equations (11) (both when a = 0, hence Sag ~ Sar = 
and not): 

r{t) ^27/5 + eacos{uj;)t+ (14) 

+ {e^p + e^pi ) cos (2^* )i + (e^c + e^a ) + 
+ e^d cos {lu* -2ujg)t + 
+ €^e cos {uj* + 2ujg)t + €^f cos {3iLj*)t + 
+ e'^r cos(2tj* + 2LUg)t + e'^scos {4:UJ*)t + 
+ e'^u cos {4:Ujg)t 
B{t) =n/2 + egcos{uj0)t + 

+ [e^h + e'^/ii) cos {u)* - uj*g)t + 

+ {e^l + eHi)cos{uj*+uj*g)t + e'mcos{2wl-ujl)t + 



-I- e^ncos(2a;* + uj*g)t + e'o cos {3ujl)t + 

+ e'*j cos(3tJ* +(jjg)t + t^w cos (uj* + 3ujg)t + 

+ e'^gcos {uj* — 3Ljg)t 

The coefficients of the trigonometric functions {a,g,c..) 
are constants which depend on rg and the initial condi- 
tions. For simplicity of notation here the initial phases 
are equal to zero. 

4. Symmetry and regions of resonance 

The Taylor expansion at constant angular momentum 
to third order near the stationary point of the effective 
potential leads to two nonlinear differential equations with 
constant coefficients: 

■£{t)+^lz{t) = a^juzity p{ty z(t)k p(ty (15) 

q^r,o.p 

The polynomials of third degree {i+ j + k + 1 = 3\/i,j,k,l) 
which constitute the nonlinear part of the equation 
depend on the symmetry of the effective potential and of 
the metric. 

The first part of the method of multiple scales consists 
in some algebraic calculations that highlight what are the 
so-called regions of resonance. Indeed the coefficients of 
the approximate solution are found to have in the denom- 
inators terms like nujr — mujg. These n and m cannot take 
any value, but they depend on the symmetry of the ef- 
fective potential, the metric and the perturbation via the 
i,j,k,l which are present in the equations : the allowed n 
and TO indentify the candidate regions of resonance. The 
symmetry restricts further these regions : it could happen 
that the constant coefficients cancel and then there would 
be no resonance at all. This calculations performed for a 
general effective potential Uef / = Uef / (r, 9) showed that 
the following resonances, usually present, cannot appear 
in the plane-symmetric case: 

2:1, 3:1, 4:1, 1:4, 1:3, 2:3 

The only possible resonances to this degree of approxi- 
mation are then: 

1:2, 1:1, 3:2 

Since in General Relativity tj,. < LUg (this is due to the 
curvature of the space), then the only possible resonance 
is the 3:2, which is indeed observed most of the times 
in black holes candidates. In the next approximation a 
new resonance appears, 4:2, and one may find many 
other possible to : n. However a property of parametric 
resonance is that the region of instability becomes thinner 
when m + n increases: the "higher" resonances are less 
probable (Bell 1957). This is why we do not observe them. 

Landau & Lifshitz ("Mechanics" 1976), when study- 
ing systems similar to ours (chapter 5 par. 28) with the 
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non-geodesic terms a = 0, use the method of Lindstedt- 
Poincarc and they conclude that these are the solution we 
arc interested in, because in a closed system, without any 
source of energy, there cannot be a spontaneus increment 
in the intensity of oscillations. In this sense we asserted 
that "no parametric resonance occurs for strictly geodesic 
motion" {Ahramowicz et al. 2003): the initial conditions 
of the perturbation of the geodesies are related to the en- 
ergy; even if the initially small two perturbations couple 
and exchange energy, their amplitudes cannot become ob- 
servable. 

5. 3:2 resonance near non-rotating compact ob- 
jects 

An approximate solution for Schwarzschild's perturbed 
geodesies was obtained and the validity of the approxima- 
tion was verified a posteriori for initial conditions of the 
order of 10~^ (which were used in the numerical integra- 
tion). 

In this solution two modes remain locked, with the two 
frequencies near to the 3 : 2 ratio: the value of these fre- 
quencies is were we would see the twin peaks in the Fourier 
power spectra of the observed signal. The correction to 
the eigenfrequencies can be written in the following form 
(higher order corrections can be neglected): 



uj* ~ Wr + t^^^n''o,a,ae(0),<5r(0)) 



(16) 



,(2) 



[ro,a,ae(0),(5r(0)] = 
1 



;;2 



+ 4:C2lUl^LOg + Cf^UJ^UjQ + Cii (-620^^,^ + 

-I- CbUj^ + 4620^9 ~ er2Lol + 4er2W,^We)] 

(0) -I- [cii(-3eo2t^^ + 8eo2t^e + 
- e22Wr'^e + 8622^^6 ) +t^r(-3co3Wr + 
-I- 12co3We - 2c;,eo2We + 2cbe^2i^g)]al{Q)} , (17) 



tj^^^ [ro,Q;,ae(0),Q!r(0)] = 

[ lOCjo + 9e30W^ + 10e2oer2t^r + 

+ 6[-ei2Wr 4- Aei2Ujlujl + 

+ eo2 (-2620^^,^ + 8e2ot^e + 
+ AcbUj^ujg) - 2e2oez2(^^uje + 



- 2cbUJgUJ^ + 8620622^^6 + 4c6e22^r^e + 
+ 2ciiU!^{eQ2-ujgez2)]al{0)} . 



(18) 

The dependence on a and 7*0 is included in the con- 
stant coefficients: one can vary the initial conditions 



(50(0) = eae{0) - z(0) ,5^(0) = ear{0) ~ p(0) ) and/or a 
^ , obtaining different values for the corrected frequencies 
and for their ratio. For small variations of these parame- 
ters, the dependence on them is almost linear. 
The analytic approximation in figure (2) was obtained by 
fixing |a| = 0.99 (far from geodesies) and ro = r3:2: each 
point on the segment corresponds to the state reached 
by the system starting from differemt initial conditions 
(r(0) G [0.23,0.21] and z{0) G [-0.26-0.19]). The same 
slope was matched numerically (in the pseudo-Newtonian 
case) by fixing the initial conditions and by varying a 
The fact that the slope is smaller than 3 : 2 depends on 
the choice of the initial conditions. This analytic result 
is important in explaining the behaviour of the solution 
more than the numerical values: the fact that the observed 
ratios are not exactly in 3 : 2 ratio is a feature of paramet- 
ric resonance and the position of the centroid frequencies 
depends on the strength and features of the non-linear 
perturbation. 

6. Conclusions 

The Taylor expansion of the relativistic geodesies to 
the third order leads to two coupled harmonic oscillators: 
the purely geodesies motion is stable, while when it is 
perturbed there are cases for which parametric resonance 
may occur. 

Using the method of multiple scales one can highlight that, 
owing to the curvature (through the way it determines the 
effective potential), the first allowed resonance between 
the radial an the vertical epicyclic frequencies is the 3:2, 
in agreement with the numerical analysis. 
Moreover the deviation of the slope from 3 : 2 is easily 
explained as a property of such a non-linear resonance 
The analysis of this toy-model reinforces the theory that 
indeed the observed pairs of QPOs may be due to para- 
metric resonance, and finally to the strong gravitational 
field alone. 
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Appendix 1. Errata corrige 

The equations studied in (Abramowicz et al. 2003) for 
the perturbed geodesies in the Pseudo-Newtonian poten- 
tial were incorrect. The correct ones are: 



ei2 = 



er2 ^ 



4 



1 



1 d d^U 
1 



.^(1-^)1 

80^ ^ ro 



eire2 



2ro ro - 
l-2ro 



1 



.ro^(ro-l)2 ' 

where the derivatives of the effective potential are eval- 
uated in the equilibrium and E = —Ut is the energy. 
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fip,z,ro,Oo) = cnzp + ci,zp + 



(Al) 



cibpzp 
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622^ 



g{p,z,ro,Oo) ^ eo2Z +e2oP'' 

esoP^ + eize2pz'^ + ei2pz'^ . 

This does not affect the numerical results: indeed the 
symmetry (then the possible resonances) is the same and 
the only difference is in the values of the initial conditions 
which lead to the resonance. 

Appendix 2. Coefficients 

In equations (11) I named the constant coefficients in a 
way practical to remember their origin : for example C21 
is a coefficient in the equation z (it has the letter c) and 
it it is the coefficient of p at the 2nd power and z at the 
first power (numbers (21) ). Every coefficient contains the 
derivatives of the effective potential and of the metric at 
equilibrium. The explicit form of each coefficient is: 



1 d d^U 



1 d^U 



(A2) 



Cb = 



C21-^'( 



3 d^U ^ J_d^dHJ_ 



24 de'' 

1 d d^U 



dr dO^ 
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C03 



12^ 



eo2 



1 . d d^U, 



ro dr de^ 



e20 = E^[ 



1 d^U 



2r2 dr^ 



1 1 d^U 
4 ^ ro' dr^^ 
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